NUMERICAL SOLUTIONS OF THE SPECTRAL PROBLEM FOR ARBITRARY 
SELF-ADJOINT EXTENSIONS OF THE ID SCHRODINGER EQUATION 
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Abstract. A numerical algorithm to solve the spectral problem for arbitrary self-adjoint ex- 
tensions of ID regular Schrodinger operators is presented. The construction of all self-adjoint 
extensions of the symmetric Schrodinger operator on a compact manifold of arbitrary dimension 
with boundary is discussed. The self— adjoint extensions of such symmetric operators are shown 
to be in one-to-one correspondence with the group of unitary operators on a Hilbert space 
of boundary data, refining in this way well-known theorems on the existence of self-adjoint 
extensions for Laplace-Beltrami operators. The corresponding self— adjoint extensions are char- 
acterized by a generalized class of boundary conditions that include the well-known Dirichlet, 
Neumann, Robin boundary conditions, etc. Only the numerical solution of ID regular cases 
are consider in this paper. They constitute however a non— trivial problem. The corresponding 
numerical algorithms are constructed and their convergence is proved. An appropriate basis of 
boundary functions must be introduced to deal with arbitrary boundary conditions as described 
by the general theory. Some significant numerical experiments are also discussed. 



In this paper we present an algorithm to solve numerically the spectral problem for all possible 
self-adjoint extensions of ID regular Schrodinger operators. The algorithm can be easily extended 
to higher dimensions, however we will restrict our attention to the ID case in this paper. As it 
will be shown later, in the ID situation, the space of self-adjoint extensions of regular Schrodinger 
operators is parametrized by a finite dimensional unitary group, while in dimension 2 and higher, 
such space is described by an infinite dimensional unitary group. Thus in the ID case, we can 
provide a complete analysis of the boundary conditions of the numerical algorithm and prove its 
convergence. We will discuss both, the singular case and the higher dimensional algorithm in 
subsequent works. 
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We will start performing the analysis of the problem in arbitrary dimension D and will restrict 
to dimension one at the moment of the development of the numerical algorithm. The study of the 
self-adjointness of Schrodinger operators has been a fundamental mathematical problem since the 
beginning of quantum mechanics and there is a large literature on the subject (see for instance 
[Re75| ■ the review [SiOO] and references therein). In spite of this, there is a continuous flow of 
results and surprises (see for instance the recent papers where some apparent paradoxical aspects of 
the spectrum of certain self-adjoint extensions of Schrodinger operator in 2D are analyzed |Be08| , 
pa09] . [Be09] 'l. 

Hence, we will consider the evolution of a quantum system on a _D-dimensional riemannian 
manifold M with boundary under the influence of a potential V which is given by the 

Schrodinger equation: 

(1.1) m—=H^S, 
with H the Hamiltonian operator of the system given by 

(1.2) H = -^A, + V, 

2m 

where stands for the Laplace-Beltrami operator on M defined by the metric f]. The second 
order differential operator H is formally self-adjoint, however in order to define a unitary evolu- 
tion of the quantum state Vt, a self-adjoint extension of it must be specified. If H denotes one 
of such self-adjoint extensions, because of Stone's theorem, a one-parameter group Ut of unitary 
operators exists such that Ut = exp{—itH/h) and the quantum evolution of a given initial state 
^0 will be uniquely determined as = Ut'^Q. Self-adjoint extensions of the operator H are 
usually chosen by fixing the boundary values of the functions where the operator H acts, typically 
Dirichlet or Neumann like boundary conditions. However they are by no means the most general 
choice of boundary conditions determining self-adjoint extensions of the operator H and a variety 
of other possibilities exist. The development of quantum information technologies makes relevant 
the discussion of more general classes of boundary conditions, i.e., of general self-adjoint exten- 
sions for the operator H, and the numerical computation of their spectrum in order to integrate 
Schrodinger's equation. 

Von Neumann's theorem (see for instance jWeSOj ) provides a characterization of all self-adjoint 
extensions of the operator H, but such approach is not satisfactory from the computational point 
of view because it involves, as a first step, the computation of the deficiency spaces of the operator 
and later on the construction of the operator on domains which are not easily described. Thus the 
current numerical algorithms used to compute the spectrum of the operator H are not well suited 
for the use of von Neumann's theorem. On the other hand, there has been significant developments 
in the analysis of self-adjoint Sturm-Liouville problems in terms of the boundary data of the space 
of functions where the operator acts (see for instance [?] and references therein). These will help 
in constructing an effective numerical algorithm to compute the spectrum of arbitrary self-adjoint 
extensions of Laplace-Beltrami operators. In this article we will take the approach developed in 
|As05| to describe the set of self-adjoint extensions of the Laplace-Beltrami operators. It was 
shown there that the set of self-adjoint extensions of the Schrodinger operator H is in one-to- 
one correspondence with the group of unitary operators on a space of boundary data for the 
problem. We review these results in Section 2 and provide proofs of the main results suited for 
the development of this paper making it, in this way, self-contained. Because of this and in order 
to make the presentation as simple as possible, we will consider that the boundary dM of the 
manifold M is smooth and that M is compact. We will also assume that the potential function V 
is regular. It is not necessary to say that the case of non-regular potentials and/or non-compact 
manifolds emerges naturally in the discussion of many examples and applications, however we will 
restraint ourselves in the analysis of this article to the simple situation mentioned above to clearly 
exhibit the main ideas and to free the discussion from technical complications. The extension of 
the results discussed here to the more general situations just mentioned will be done elsewhere. 



^Regularity conditions for the manifold and its boundary will be discussed later on. 
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The numerical algorithms we will discuss here are based on the finite elements method commonly 
used to solve elliptic problems in riemannian manifolds. Because we are going to construct the 
appropriate finite elements to handle general boundary data, the algorithms can be developed 
provided that the boundary manifold dM can be characterized nicely. If the dimension D oi M 
is 1, 2 or 3, the boundary dAI will be respectively 0, 1 or 2 dimensional. In all these cases we 
have a complete characterization of such manifolds. Thus in the ID case, the boundary manifold 
will consists on a finite family of points. In the 2D situation, the boundary dM will consists on 
a finite collection of oriented circles and in the 3D case, the boundary will consists on a finite 
family of compact oriented Riemannian surfaces. Notice that if we were considering non-compact 
manifolds, the number of connected components of dM could be infinite. Moreover, in such a 
case singularities could occur. For instance, some connected components of the boundary could 
collapse to points in the 2D or 3D situations, etc. In this paper we will discuss in detail just the 
ID case because the boundary of the manifold has the simplest structure. Similar analysis can be 
conducted in the 2D and 3D situation by using the appropriate discretizations of the boundaries, 
but we will not dwell into such analysis here. 

The paper is organized as follows. Section 2 will be devoted to prove the general theorems that 
will allow for a convenient description of the self-adjoint extensions of the Schrodinger operator 
H in terms of boundary data in arbitrary dimensions. In section 3 we describe the particular form 
that such theorems take in ID and we prepare the setting for the development of the corresponding 
computational tools. The main results of the finite element method for the eigenvalue problem 
in the interval, for general self-adjoint extensions, are discussed in section 4. In section 5 some 
numerical experiments displaying some relevant features of the algorithms are shown. 



2. Self-adjoint extensions of Schrodinger operators 



As it was stated in the introduction we will restrict our attention to the case of Schrodinger 
operators on compact manifolds with smooth boundary and regular potentials. The Schrodinger 
operator for a particle moving on a smooth manifold M with boundary dM and riemannian metric 
7] is given by the Hamiltonian operator H defined in eq. (1.2) that, in local coordinates x% takes 
the form: 

1 d 



2m 



with the metric tensor rj given by rj 
second order differential operator 

(2.1) = 

is formally self-adjoint in the sense that 



mv'''— 

7y| dx^ dx^ 

rijk{x)dx^dx^ , \ri\ = 



I detrijk{x)\ and r/*-'??jfe 



The 



1 



d 

dx^ 



■,k 9 
'^l^ dx-^ 



(2.2) 



M 



A^<I> vol„ 



M 



M 



(A^4')$vol^, 



for $ any smooth complex valued functions with compact support contained in M\dM. In 



the previous formula, eq. (2.2), vol^ denotes the riemannian volume form on M defined by rj, 
^y\r]\dx^ A • • • A dx^, and (d^I', o?$)^(2.) is the inner product among covectors at a; e M defined by 
the expression rj^'^ {x)d^ / dx^ d^ / dx'^ . In fact, the differential expression (2.1 ) defines a symmetric 



operator on the space L'^{M) of square integrable functions on M with respect to the measure 
defined by the volume form vol^, with dense domain C^{M\dM) the space of smooth functions 
with compact support on the interior of M. The operator A.,, is symmetric but not self-adjoint on 
L'^{M) because the adjoint operator A^ has a dense domain ~ H^{M), the space of functions 
of Sobolev class 2 on M, i.e., the space of functions on L'^{M) possessing first and second weak 
derivatives that are square integrable. Clearly {M\dM) C H'^{AI). The operator A,, is closable 
on L'^{M) and its closed extension, that we will keep denoting by A^, has domain Dq :— Hq{M), 
the space of functions of Sobolev class 2 on M vanishing at the boundary. Self-adjoint extensions of 
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the operator are given by operators A^i with domain D such that Dq C C -DJ, A^j |_Do== A,, 
and Ax) — A|j. Notice that in that case A/5 = AjJ/j. 

Von Neumann's theorem estabhshes (see for instance |We80| . Thm. 8.12) that there is a one- 
to-one correspondence between self-adjoint extensions A^ of the Laplace-Beltrami operator A^ 
and unitary operators K : Af+ — > A/L, where the deficiency spaces A/± are defined as: 

AA± = e H^{M) I A],* = ±i^}. 

In particular, given the unitary operator K, the domain D of the operator A^i is given by D = 
Dq (B {I + K)M+ , and the extended operator A £> takes the explicit form: 

Az,(*o © (/ + K)^+) = A^^-o ® til - K)^+, 

for all ^'o e Do and ^+ G 7V+- 

Unfortunately, as it was stated in the introduction, von Neumann's theorem is not always 
suitable for the explicit construction of general self-adjoint extensions of the Laplace-Beltrami 
operator because we need to determine first the deficiency spaces Af±. We can take however a 
different route inspired in the classical treatment of formally self-adjoint differential operators. 



If we rewrite eq. ( [2^ for functions $ in ZJ'f = 'H'^{M) instead of C^{M\dM), a simple 
computation shows: 

(2.3) / ^'A^$vol„= / {K^)^Yo\+ j U^-i^^)Yo\an. 

JM J M JdM ^ ' 

where '^l) — \qm t ^ ~ ^ Isa/j f^nd the normal derivative ^ is defined as: 

^vola^ = *d<I> \qm 

where -k is the Hodge operator defined by the metric r\ij and volg^ is the riemannian volume defined 
on the closed boundary manifold dM by the restriction dr\ of the riemannian metric 77 to it. Less 
intrinsically, but more explicitly, we have 'P = ^ \dM= d^{i') where ly is the exterior normal 
vector to dAI. We thus obtain the Lagrange boundary form E for the Laplace-Beltrami operator: 

(2.4) j:{{tp,ip),{(p,(p)) = {ip,(p)mdM) - {i^,'p)LHdM)- 

In what follows, if there is no risk of confusion, we will omit the subscript L'^{dM) that denotes 
the inner product on the boundary manifold dAI with respect to the measure defined by the 
volume form volg,;, hence we will simply write {ip,ip) = J^j^ rpipvoldri- The Lagrange boundary 
bilinear form S defines a continuous bilinear fom on the Hilbert space L^{dM) © L'^(dM), 

(2.5) E((7/^i,7A2),(¥'i,¥'2)) - (^1,(^2) - {^2,m), V(^i,, 1^2), (^1,^2) e L\dM)(SL\dM). 

If we denote by 7: 'H'^{M) L^{dM) ® L'^{dM) the trace map given by 7(^') = {ip,^^), 7, it is 
continuous and induces a homeomorphism from ?^^(M)/kcr7 into W^^^'^{dM) W^^'^''^{dM) C 
L^{dM) ® L'^{dM) (see for instance |Ad75j . thm. 7.20). The previous observations provide a 
simple characterization of self-adjoint extensions of the operator A^. In fact it is easy to check 
that: 

Theorem 2.1. There is a one-to-one correspondence between self-adjoint extensions A]j of the 
Laplace-Beltrami operator A^ and non-trivial maximal closed isotropic subspaces W of Y] con- 
tained in W^^'^''^{dM) W^^'^''^{dM). The correspondence being explicitely given by D l{D). 

Proof. Let D C 'H^(M) be the domain of a self-adjoint extension A^) of the operator A,,. Consider 
the subspace W := 7(£i) C W^''^-^{dM) ® W^/^^'^{dM) C L'^{dM) ® L^{dM) consisting on the 
set of pairs of functions (ip, ip) that are respectively the restriction to dM of a function <& G _D 
and its normal derivative. Notice that the subspace D is closed in 1-1^ {M). Hence, because 7 is an 
homeomorphism from n'^{M)/keT-f to W^I'^J^M) ® W'^/'^^'^{dM), 7(D) is a closed subspace of 



W^/'^^'^{dM)(BW'^/'^^'^{dM). Because of eq. ([2^, it is clear that E and finally the subspace 
W is maximally isotropic in iy'^/^'^(9M) ® W^^'^(9M), because if this were not the case, there 
will be a closed isotropic subspace W C W^/'^''^{dM)®W'^/'^^'^{dM) containing W. Then 7-1 (W) 
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will define a domain D' containing D such that the operator A,, would be symmetric on it, in 
contradiction with the self-adjointness assumption. 

Let W C L'^{dM) ® L'^{dM) be a maximal closed E-isotropic subspace in Ty^/^^^(9M) © 
W^i/2'2(aM). Then consider the subspace Dw ■= I'^iW) =C 'H'^{M) of functions ^' such that 
(tpjip) G W. It is clear that for any pair of functions $ on Dw, because W is isotropic with 
respect to S, then eq. (2.3 1 gives (^', A^$) = (A^^f, $) and the operator A^ is symmetric in Dyy- 
Moreover, because of the maximality of M^in W^^'^''^{dM) © W^/^'^(9M) it is easy to see that 
Dw — d\^, hence it is self-adjoint. □ 

We can object that the previous characterization of self-adjoint extensions of the Laplace- 
Beltrami operator in terms of closed maximal isotropic subspaces of S contained W^^'^''^{dM) ® 
Ty^/^'^(9Af) is rather obscure and not easy to describe explicitly. An important observation in 
this sense is that the linear transformation C: L'^{dM) ® L'^{dM) L^{dM) © L'^{dM), defined 
by C{ip, (f) = + -^{ip — iip)), transforms maximally isotropic closed subspaces of S into 

graphs of unitary operators of L'^{dM). 

Theorem 2.2. The map C provides a one-to-one correspondence between maximally isotropic 
closed subspaces of the Lagrange bilinear boundary form S in L^{dM) © L'^{dM) and graphs 
Gv ^ {i^p+,ViJ+) I ^p+ G L^{dM)} of unitary operators V : L^{dM) L^{dM). 

Proof. Notice first that the map C is a unitary operator on the Hilbert space (dM) © (DM) 
with the inner product (((f/'i, V'2), iipi,ip2))) = + {'^2,^2)- 

Consider now the transformed bilinear form E on L^{dM) © L^{dM) defined by 

E(c(^, ^), V')) = m^p, v), (V', V')). 

Thus, using the notation {p± — ^{^p ± iip), we have: 

S((v5+,<^-), (^+,■0-)) = [(<^+,V'+> - (■^-,-0-)] • 

Hence, if is a maximally isotropic closed subspace for S, then W = C{W) will be a max- 
imally isotropic closed subspace for E. Then it is easy to show that W defines the graph of 
a linear operator. We first realize that ({0} x L^) n = because if {0,ip-) lays in W, then 
5i((0,V'-), (0,-0-)) = must vanish. Hence if (-0+ , V"- ) ; (V'-i-i'^A-) are in W', then {0,tp_—'tp_) 

is in W and xjj^ ~ Then we define a linear operator V:V ^ l? hy U{ip+) — ijj- with 

(V'-i-, V'-) G W and V the closed subspace of vectors -0+ such that there exists {4'+, V'-) G W . Sim- 
ilarly we can construct another operator y : V — ^ by observing that (L^ x {0}) C\W = 0. Then 
it is easy to show that V is an isometry from V to V, is an isometry from V to V and they are 
inverse of each other. Then because of the maximality of W , we conclude that V = V = . □ 

Hence a convenient way of constructing self-adjoint extensions of the Laplace operator will be 
provided by unitary operators on L'^(dM) such that the preimage under C will be contained and 
closed in W^/^''^{dM) © W^/'^''^{dM). We will develop this programme fully in the forthcoming 
section in the ID case. 

We will end this discussion by realizing that the operator multiplication by a regular function 
is essentially self-adjoint and its unique self-adjoint extension has domain L^(M). Hence, the 
self-adjoint extensions of H coincide with the self-adjoint extensions of A^. 

We can summarize the preceding analysis by stating that under the conditions above, the 
domain Z) of a self-adjoint extension of the Schrodinger operator H is defined as a closed subspace 
of functions on TH?{M) satisfying: 

(2.6) ilj-ii! = V{ip + ii!) 

for a given unitary operator V : L^{dM) — > L^{dM). In what follows we will denote respectively 
by Hy or Hjj the self-adjoint extension determined by the unitary operator V or the self-adjoint 
extension whose domain is D. Also according to Thms. 2.1 and 2.2 we denote such domain as 
Dv- Notice for instance that V = I corresponds to Neumann's boundary conditions and V = —I 
determines Dirichlet's boundary conditions. 
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The formula above, eq. (2.6), provides a powerful and effective computational tool to deal with 
general self-adjoint extensions of Schrodinger operators. It was introduced in a slightly different 
context by Asorey et al |As05) and will be used extensively in the rest of this paper. 



3. Self-adjoint extensions of regular ID Schrodinger operators 

3.1. The unitary group of self— adjoint extensions in ID. Instead of continuing the dis- 
cussion about self-adjoint extensions of Schrodinger operators in arbitrary dimensions we will 
concentrate our attention in ID. Then we will be able to provide an elegant algorithm to solve the 
spectral problem for each self-adjoint extension. 

Notice first that a compact ID manifold M consists on a finite number of closed intervals /q,, 
a = 1, . . . , rt. Each interval will have the form — [oq, ha] C M and the boundary of the manifold 
M = \XL^i[(ia,ha\ is given by the family of points {ai, 5i, . . . , a„, 6„}. Functions ^ on M are 
determined by vectors (^'i, . . . , ^'„) of complex valued functions ^"0, : — C. 

A riemannian metric 77 on M is given by specifying a riemannian metric rja on each interval la, 
this is, by a positive smooth function ria(x) > on the interval la, i-C, rj |/^= ria{x)dx'^. Then 
the inner product on takes the form (^I'q,, = {x)^a{x)^aix)y/r]a{x)dx and the Hilbert 
space of square integrable functions on M is given by L'^{M) = ®a=i L^^ajVa)- The Hilbert 
space L^{dM) at the boundary reduces to vectors in C^", as well as the subspaces W'^^^'^ldM) 
and W^^^'^{dM), determined by the values of ^I^ at the points aa, ba with the standard inner 
product: 

ip = (*i(ai),*i(foi),... ,*„(a„),^'„(fe„)). 
Similarly we will denote by ip the vector containing the normal derivatives of at the boundary: 









I dx 


ai 


dx 



bi 



dx 



d^r 



dx 



Because of Thm. |2.2| an arbitrary self-adjoint extension of the Schrodinger operator H defined 
by the riemannian metric 77 and a regular potential function V is defined by a unitary operator 
V : C^" — ^ C^". Its domain consists of those functions whose boundary values ip, tp satisfy equation 



(2.6) above. This equation becomes a finite dimensional linear system for the components of the 
vectors ip and ip. Hence the space of self-adjoint extensions is in one-to-one correspondence with 
the unitary group U{2n) and has dimension 4n^. 

It will be convenient for further purposes to organize the boundary data vectors ip and ip in 
a different way. Thus, we denote by ipi G C" (respec. ipr) the column vector whose components 
ipi{a), a = 1, . . . ,n, are the values of ^' at the left endpoints this is ipi{a) — ^^aio-a) (respec. 
ipr{oi) — ^q(&q) are the values of 5* at the right endpoints). Similarly we will denote by ipi{a) = 
— Iqc ^^'^ ^r(a) — '^%pr l&Q, 01 — l,...,n. Hence, the domain of the self-adjoint extension 
defined by the unitary operator V will be written accordingly as: 



(3.1) 



ipi-iipi = U^^{ipi+iipi) + U^'^(ipr + iipr) 

Ipr-ilpr = U'^\lpl+ilpl) + U'^'^{lpr+ilpr) 



where the unitary matrix U is obtained from V G U{2n) by the permutation a 
2 3 n + 3 • • • n 2n) and has the block structure: 



(1 n + l 2 n- 



(3.2) 



U 





JJ12 - 


JJ21 


JJ22 



Thus in what follows we will use the notation for the boundary data: 



iP 



ipi 

Ipr 



iP 



Ipl 
Ipr 



and the boundary equation (2.6) is written now as: 
(3.3) 



iP-iiP = U{ip + iip), UeUi2n). 
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3.2. The spectral function. Once we have determined a self-adjoint extension Hjj of the 
Schrodinger operator iJ, we can determine the unitary evolution of the system by computing 
the flow Ut = exp{—itHu/h). It is well-known that the Dirichlet extension of the Laplace- 
Beltrami operator has a pure discrete spectrum because of the compactness of the manifold and 
the ellipticity of the operator, hence all self-adjoint extensions have a pure discrete spectrum (see 
[We80| . Thm. 8.18). Then the spectral theorem for the self-adjoint operator Hu states: 



Hr. 



k=l 



where Pk is the orthogonal projector onto the finite-dimensional eigenvector space Vk correspond- 
ing to the eigenvalue Afc. The unitary flow Ut is given by: 



k=l 

Hence all that remains to be done is to solve the eigenvalue problem: 

(3.4) iJ,7* = A*, 

for the Schrodinger operator Hu . We devote the rest of this paper to provide an efficient numerical 
algorithm to solve eq. (3.4). However, before embarking into that, it is remarkable that the ID 
situation allows for an explicit formula that determines the spectrum of Hjj and that can be used 
also to compute numerically the eigenvalues of Hu. We shall derive such formula first. 

On each subinterval = [acfoa] the differential operator Ha = H\j^ takes the form of a 
Sturm-Liouville operator 

W a dx ax 

with smooth coefficients Wa — ^ (now and in what follows we are taking the physical 

constants h and m equal to 1), Pa{x) = :y|=j hence the second order differential equation 

(3.5) Ha-^a - A*, 

has a two-dimensional linear space of solutions for each A. We shall denote a basis of solutions of 
such space as ^'J^, tr = 1,2. Notice that depends differentiably on A. Hence a generic solution 
of eq. (3.5) takes the form: 

Aa,i'fl + Aa,2'fl. 



Now it is clear that 



Hence: 



Ma) = *a(aa) = Aa^i^plia) + Aa^2i'lia). 



where Arr, a 



1,2, denotes the column vector 

" Ai^a 

A,= : 

An,c 

and o denotes the Hadamard product of two vectors jHo94j . i.e., {X oY)a = XaYg where X, y e 
C". We obtain similar expressions for ipr, tpi and 4'r- With this notation eqs. (3.1) become: 



(3.6) (V/ 



''0A2 



ii)])oA^ 

ii>\)oA^ 



u 



m) ° A2 

#?) o A2 



It will be convenient to use the compact notation ip^ 



i± 



+ iijj})oAi 
+ i-^l)oAi^U'^^{iPl 

± iipi^, a = 1,2, and similarly for ip', 



ii'r) ° A2 



r±- 
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If T is a n X n matrix and X, Y arbitrary n x 1 vectors, we will define T o X as the unique 
matrix such that (T o X)Y = T{X o Y). The rows of the matrix T o X are o X or alternatively, 
the columns of T o X are given by Xj (no summation on j). It can be proved easily that 

(3.7) ToX ^To{X(g)l), 

where 1 is the vector whose components are all ones (i.e., the identity with respect to the Hadamard 
product o) and the Hadamard product of matrices in the r.h.s. of eq. (3.7 1 is the trivial compo- 
nentwise product of matrices. Using these results eqs. (3.6) become: 

{In o ipl- 







c/ii o v-iV - u'^ o + (/„ o v-; 

(J„oV'l_-i721oV/+-t/22o^l+)Ai + (/„0 7A^^_-[/^^o^f+-C/^^o^^\)A2 = 



^11 

r21 



t22 



Thus the previous equations define a linear system for the 2n unknowns Ai and A2. They will 
have a non trivial solution if and only if the determinant of the 2n x 2n matrix of coefficients 
M{U,X) below vanish: 



M{U, A) = 



/„ O V-^^ - U^^ O _ ^22 „ ^2 



The fundamental matrix M(U, A) can be written in a more inspiring form using another oper- 
ation naturally induced by the Hadamard and the usual product of matrices. Thus, consider the 



2n X 2n matrix U with the block structure of eq. (|3.2|) and the 2n x 2 matrices: 

[^i I i^l] 



then we define 





yl2 ■ 


C/2i 


C/22 _ 









1 






°^]± 


+ [/12 






o ijf^ - 


|-[/12 














+ C/2^ 




^2i 






o V4 . 



and similarly 



hn [V'i I V'd 



Finally we conclude that the condition for the existence of coefficients Ai and A2 such that the 
solutions to the eigenvalue equation lie in the domain of the self-adjoint extension defined by U is 
given by the vanishing of the spectral function Ay (A) = detM{U, A), or written with the notation 
introduced so far: 



(3.8) 



Au{X) = det(/2„ [V'i I V-] - C/ [V'i I ^+]) - 0. 



The zeros of the spectral function A provide the eigenvalues A of the spectral problem eq. (3.4 1. 

In the particular case n — 1, the previous equation becomes greatly simplified, the Hadamard 
product becomes the usual scalar product and the Hadamard-matrix product is the usual product 
of matrices. After some simple manipulations, the spectral function Au(X) becomes: 

(3.9) Au{X) = W{l,r,-,-) + U^^W{r,l,-,+) + U^^W{r,l,+,-) + U^^W{r,r,-,+) 

+U^^W{1, 1, +, -) + dct [/ • W{1, r, +, +) 

where we have used the notation: 

^1+ 



Wil,l, 



^1 



^1+ 



V>2 

4>l 



1+ 
2 



etc. 



If we parametrize the unitary matrix U £ U{2) as: 



U 



a_ P 

~P a 



then the spectral function becomes: 
(3.10) Au{X) = W{l,r,- 

-pw{i.i,+,-) 



) + aW{r, /,-,+) + aW{r, I, +, -) + /3T^(r, r, -, +) 



Wil,r, 
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In particular if we consider a s ingle interval [0, 27r] with trivial riemannian metric, the fundamental 



solutions to the equation eq. (3.51 have the form 'i'^ — e 

W{l,r,-,- 
W{l,l,+ 
W{r, r, — 



/2\ 



^ and 



^. Then we have: 



= -2i{l + 2A) sin(27rV2A) - 4V2A cos(27rV2A), 

= 4v^, 

= 4%/2A, 

W{r,l,-,+) ^ 2i(l - 2A)sin(27rV2A), 

W{r,l,+,-) = 2i(l - 2A)sin(27r%/2A), 

W{l,r,+,+) = -2i(l + 2A)sin(27r\/2A) +4\/2Acos(27r\/2A), 

and finally we obtain the spectral function for this simple situation: 

A,7(A) = -2i(l + 2A)sin(27rA/2A)-4\/2Acos(27r\/2A)+4iRe(a)(l-2A)sin(27r\/2A) 
+8Im(/3)\/2A + e'^[-2i(l + 2A) sin(27r\/2A) + 4\/2A cos(27r\/2A)]. 



4. Finite element method for the eigenvalue problem in the interval 

4.1. Weak formulation of the spectral problem. Summarizing the discussion above, we con- 
clude that solving the evolution problem for a given self-adjoint extension of the Schrodinger 
operator on the ID manifold M — Y[^=i['^a,ba], amounts to solve the corresponding spectral 
problem determined by a unitary matrix U explicitely given by the differential system: 



(4.1) 



- A,,$ + y$ = A$, $ e -H^ (M) 

If — i(p = U {f + i(p) , (p — ^\ 



IdAP 



dv IdM' 



U e U{2n). 



As before we suppose that the potential V is regular. We will develop in the following paragraphs 
a numerical algorithm to solve such system. 

From now on we will keep the dot notation for the outwards normal derivative at the boundary 



d<S> 
du 



dM 



and use primes to denote the standard derivative <i>' = 4^ . We will denote by a "half" 



of the Lagrange boundary form (2.4 1, this is: 



Because of the particularly simple form that the Lagrange boundary form takes in ID, it will be 
also convenient to introduce the following notation: 

n 

[^maM = <J{{^, iv, ^)) = J2 ^(^a)<i>'(&a) - *(aa)<i>'(a„). 



In order to design an algorithm for solving this problem, it is convenient to look for weak 



solutions of (4.1 ). Taking the inner product of (4.1 ) with a vector ^E" on the dense domain C°°{M) 
of smooth functions on M, and integrating by parts we will obtain: 



(4.2) 



Proposition 4.1. The function $ is a weak solution of the eigenvalue problem, i.e. it solves (4.2 I 



if and only if it solves de Schrodinger spectral problem (4.1). Moreover, if "if and $ satisfy the 
boundary conditions 



(4.3) 

the bilinear form 
is hermitian. 



If — i(p — U{ip + iip), 

g(^',$) = (^'',$')-[**']9A/4 
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Proof. First we prove the second assertion. The forms and (^|y<f>) are triviahy hermitian, 

so it is only necessary to prove that 



dM- 



However this is true as long as $ and 4" satisfy tp — iip = U{if + iip), because that is the condition 
for the maximally isotropy of subspaces for the Lagrange boundary form 

S]((V',V'),((/?,(/3)) = (V',</5)l2(9M) - {■^,V)L^{dM) = [*$']9M - W^dM- 



If <I> solves the spectral problem (4.1 1, just by taking the inner product in L\M) with e C°°{M) 



and integrating by parts once, one recovers the weak form of the problem (4.2). On the other 



hand, starting with (4.2), integrating by parts using (2.3) one has 
(4.4) 



(-A^ + y - A)$) 



0. 



But the last equation holds for all \I/ G C°°{M), i.e. in a dense subset of L'^{M), so (4.1) is 
satisfied. □ 



4.2. Finite elements for general self adjoint boundary conditions. To solve (4.2) we use 



a Ritz-Galerkin approximation to the eigenvalue problem. In order to do this we will construct a 
family of finite dimensional subspaces of functions of L'^{M) satisfying the boundary conditions 
(4.3). Such finite dimensional subspaces will be constructed using finite elements. The finite 



element model {K, V,N) that we use is given hy K = [0, 1] the unit interval, V the space of linear 
polynomials on K and N the vertex set {0, 1}. 

The domain of our problem is the manifold M which consists on the union of the intervals 
la = [o-on ba\, ct — 1, ... ,71. For cach N we will construct a non-degenerate subdivision as 
follows. Let Ta be the integer defined as = [LaN/L] + 1, where [a;] denotes the integer part of 
X. We will assume that each > 2, and N > 2n. Let us denote by r the multi-index (ri, . . . , r„). 
Then |r| = ri + • • • + r„ satisfies: 

N < \r\ < N + n. 



Now we will subdivide each interval la into + 1 subintervals of length: 

La 



h„ = 



(It could be possible to use a set of independent steps ha , one for each interval, however this would 
create some difhculties later on that we prevent in this way) . Each subinterval la contains + 2 
nodes that will be denoted as: 



„(") 



~\~ h(j^k^ k — 0, . . . J 7*0; ~\~ X. 




Figure 1. Bulk function at the i-th node 
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4.2.1. Bulk functions. Consider now the family J-r of |r|— 2n piecewise linear functions {fk'\x)}k'L2^ 
that are zero at all nodes except at the fc-th node of the interval , where it has value 1 (see fig. 
[T]), or more explicitly (a = 1, . . . , n, 2 < fc < Tq, — 1): 

S X ~ + S^ai < s < 1, 

ff'^^)^\ 1-s = < s < 1, 

otherwise. 

Notice that these functions are differentiable on each subinterval. All these functions satisfy 



trivially the boundary conditions (2.6 1 because both the functions and their normal derivatives 



vanish at the endpoints of each interval. We will call these functions bulk functions. 

4.2.2. Boundary functions. We will add to the set of bulk functions a family of functions that will 
implement non-trivially the boundary conditions determining the self-adjoint extension. These 
functions will be called boundary functions and the collection of all of them will be denoted by S,-. 
Contrary to bulk functions, boundary functions need to be "delocalized" so that they can fulfill 
any possible self-adjoint extension's boundary condition. 




Figure 2. Boundary function jS^'" 



Because the endpoints Xq 

J") 



(a) 



(a) 



ba of the intervals and the adjacent nodes, 



and XrJ, are going to play a prominent role in what follows, we will introduce some notation that 
will take care of them. We will consider the index I — 1, ... , 2n; I odd is given by Z = 2q; — 1, and / 
even by I — 2a. Now for each vector w = (wi) e C^" consider the following functions (see fig. [2|: 

f2Q-l + s(w2a_l - W2a-l) X = X^^"^ + sha, 0<S<1, 



W2a-l(l-s) X 
SW2a X 
W2a + S{v2a " W2a) X 



1 < a < n. 



XrJ + Sh 

a 1 



< X < X 



(a) 



Each function on the previous family is determined (apart from the vector w) by the vector 
V = {vi) £ C^" that collects the values of /S'™) at the endpoints of the subintervals. We select 
now a family of vectors w sufficient to span a linear space that contains a solution to the 



spectral problem with boundary conditions (2.6). Thus we denote by w^'^ the vector such that 
Wj*'' = Sii, i — l,...,2n, i.e., the 2n vectors w^'^ are just the standard basis for C^". The 
corresponding functions /J^™^ will now be denoted simply by Z?*-*-* and they have the form shown 
in figure [sj Notice that each boundary function I3^^\x) is completely characterized by the unique 
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non extremal node where it doesn't vanish and the values at the endpointsj^ We denote by v- 
I = 1, . . . , 2n, the boundary values of the functions /J*^*' above. 




Figure 3. Two different boundary functions of S 



N 



4.2.3. The boundary matrix. The 2n extremal values vi^'^ of the boundary functions /S*-*-* are unde- 



fined, but we are going to show that the 2n conditions (4.31 imposed on the boundary functions 



constitute for them a determinate system of linear equations. Because the boundary functions are 
constructed to be piecewise linear, the normal derivatives of these functions at the boundary can 
be expressed as follows. For the left boundaries of the intervals, i.e., at the points a^, we have: 



(4.5) 



dn 



I l^2a-l '^2q-iJj 



and for the right boundaries we have: 



(4.6) 



dn 



Thus the vector containing the normal derivatives of the function /S*-*-*, consistently denoted by 
is given by: 



(4.7) 131' 
where we use again the consistent notation /i/ = /i^, if Z = 2a — 1, or Z = 2a. For each boundary 



function /3^'^' , the boundary conditions (4.3) read simply as the system of 2n equations on the 
components of the vector w^'^: 

or, componentwise: 

2n 

Collecting coefficients we get 



/ = 1, . . . ,2n. 



2n 

E 



2n 



Si J - Uij 



■^If i = 2q! — 1, it is the node and if i = 2a, it corresponds to the node x^,'^ . 
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and substituting in the expressions it;|*' = Sn we finahy obtain 



2n 

(4-8) E 



hi hi 



This last equation can be written as a matrix equation 
(4.9) FV = C, 

with V a,2nx2n matrix whose entries are given by Vji = , i, j — 1, . . . , 2n. The i-th column of 
V contains the boundary values of the boundary function (3^-^\ The 2n x 2n matrix F with entries 

will be called the boundary matrix of the subdivision of the domain M determined by the integer 
N, and 



are the inhomogeneous terms of the linear systems (4.9), one for each boundary function. Using 
a compact notation we get: 

F = diag(l - i/h) - ;7diag(l + i/h), C + C/)diag(l/h), 

where 1/h denotes the vector whose components are 1/hi. Notice that F depends just on U and 
the integer N defining the Ritz-Galerkin approximation to the weak spectral problem. 



4.3. Conditioning of the boundary matrix. We will study the behavior of the system (4.9 1 
under perturbations, in other words, we compute the condition number of the boundary matrix 
F and show that it is low enough to assure the accuracy of the numerical determination of our 
family of boundary functions /3^'^. The relative condition number we want to compute is 

>,{F)^\\F\\-\\F-% 

In our case the boundary matrix F can be expressed as 

F = D - UD = {I - UDD-^)D, 

with Djj(h) — Dij = (1 + -j^Sij). Notice that the product UDD^^ is a unitary matrix which we 
will denote as Uo{h) or simply Ua if we don't want to emphasize the h dependence of Uq- Thus, 
F={I- Uo)D and 

||F|| = ||(/-{/o)5||<||/-C/o||-||5||<2||i?||. 

On the other hand 

\\F-'\\ = \\D-\l-Uo)-'\\ < \\D-'\\ ■ \\iI-U„r'\\ = ... 

minAespcc((7o){|l ^ M\ 

and thus we obtain: 

k{F) < n{D)- — . 

As D is a diagonal matrix its condition number is given by 



+ 1 h 



with /imax (^min) the biggcst (smallest) step of the discretization determined by N . We get finally, 
(4.10) «(F)<|^-^, 

with A the closest element of the spectrum of to 1. Of course, because Uq is unitary, it may 
happen that 1 is in its spectrum so that the condition number is not bounded (and the system 
(4.9) will be incompatible). Because the matrix Uq depends on h, its eigenvalues will depend on 
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h too. We want to study the dependence of the closest eigenvalue to 1, or 1 for that matter, with 
respect to perturbations of the vector h. 

Lemma 4.2. Suppose is an eigenvector with eigenvalue 1 of Uq, then if the perturbed matrix 
U — Uq + SU , for \\SU\\ small enough, is such that 1 € <^{U), then X^SUXq = Q to first order in 

su, (x^ = xi). 

Proof. Clearly, if 1 G o'(U) and SU is small enough, it exists a vector X = Xq + 6X, with 
\\SX\\ < C\\SU\\, such that UX = IX. Then, 

{Uo + SU){Xo + 6X) =Xo + SX 
^ UqSX + 6UXa + SUSX = SX. 

Because UqXq — Xq and Uq unitary, X^Uo ~ , then multiplying on the left by Xq and 
keeping only first order terms, we get the desired condition: 

X^SUXo = 0. 

□ 

Because of the previous lemma, ii U ^ Uq + SU is a unitary perturbation of Uq such that 
X^SUXq 0, for any eigenvector Xq with eigenvalue 1, then 1 ^ cr{U). Now if we consider a 
unitary perturbation U of Uq such that 1 ^ cr(?7) we want to estimate how far away 1 is from 
the spectrum of U. Consider the eigenvalue equation for the perturbed matrix. The perturbed 
eigenvalue X = 1 + SX will satisfy: 

(4.11) {Uq + SU){Xq + SX) = (1 + 5X){Xq + SX) 

Multiplying on the left by X^ and solving for \6X\ it follows that, 

^ iXj^SUXQ+Xj^SUSXl ^ \X,"SUXq\ - ixj^susxi 

- \i+x,"sx\ - rnra ' 

for ||^J7|| small enough. Considering the particular form of the matrix Uq = UDD~^ we have 
that SU = US{DD-^) and therefore \\SU\\ ^ \\5{DD-'^)\\. Having into account that S{DD-^) 
is a diagonal matrix, in fact S{DD^^)kk — (^ii^-\)2 Shk, its singular values are the modulus of its 
eigenvalues, and we have the following proposition. 

Proposition 4.3. Given the matrix Uq = UDD~^ with eigenvalue 1, then 1 is not an eigenvalue 
of the perturbed matrix U + SU = Uq{DD~^ + S{DD~^)) with \ \S{DD~^)\\ small enough, and the 
perturbation SX of such eigenvalue satisfies the following bound, 

a„,,^{6{DD-')) - Cal,,{S{DD-^)) 



\5X\> 



Ca,^US{DD-^)) 



Proof. Because of Lemma 4.2 it is sufficient to show that X^ 5{DD ^))Xq ^ 0. But this is an easy 
consequence of the fact that \ Y.IZi{\Xq^S{DD-^))u\ > 2(5/imi„. Furthermore \\SX\\ < C\\SU\\ 



and the bound follows from equation (4.12). □ 



If we neglect terms \hf\ <^ 1 and \Shi\ ^ \hi\ we finally get the desired bound for the condition 
number: 

(4.13) ,(F)<^i-, 



where Sh — Tahi{Sha}. Then, if for a given N we obtain a boundary matrix F which is bad 
conditioned, it suffices to change the size of the discretization, i.e., to increase N, to improve the 



condition number. Of course, if N is already quite large, then the bound (4.13) could be useless. 
For typical values h « 10^^ ^ 10^"^, it can be taken Sh « 10^^ ^ 10^^, to provide condition 
numbers k{F) « 10^ ~ 10^. 
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4.4. The spectral pencil. For any N > 2n we will define the finite-dimensional Ritz-Galerkin 
approximation space as the linear span of the of bulk and boundary functions, i.e., = 



span{/,",/3' |« = 1, 



, ri, fc = 2, . . . , Tq, — 1, / = 1, . . . 2n}. All functions /j^""* and /3^'^ are linearly 

It is convenient to organize the 



(2n-l) An) 



independent, thus the dimension of 5*^ will be \r\ = ri 
basis above as follows: 

Using this ordering, we will rename the elements of this basis as /a, with a = 1 
arbitrary element <&Ar G will be written as: 



and an 



We consider now the approximate eigenvalue problem: 

(4.14) («'^, - [*Ar$^]aA/ + (*7V, V<i>N) - AAr(*jv, *jv) 



0, V^-jv e S 



N 



Introducing the expansion above in the approximate problem (4.14) we get: 
kl 

(4.15) ^^(^)) + [/a(^)/^(^)]9A/ + {fa{x), V h{x)) - \{fa{x), h{x)) 
a,b 

As ( 4.15[ ) holds for every "^n & , this equation is equivalent to the matrix pencil 

(4.16) A* = AS*, 

where ^afc = {f^ix)j;,{x)) + [fa{x)f;,{x)]dM + {faix),Vfbix)) and Bab = {fa{x)Jb{x)). Notice that 
A and B are hermitian matrices, what improves the algorithms for finding the pencil's solutions. 
In fact, when solving numerically the spectral pencil (4.16), it is relevant to keep the pencil's 
hermitian character. Notice that the boundary functions /?^'' {x) satisfy 

(4.17) [;3(')/3(")']aA/ - [/3(")^('^']9M = 0, 

because their boundary values are elements of a maximally isotropic subspace of the Lagrange 
boundary form. Using (4.7) and the definition of the boundary values of the boundary functions 
in terms of V we have that: 



OM 



k=l 



2n 

fe=i ' 



-V, 



Till ■ 



This identity together with equation (4.17) leads to 

1 



(4.18) 



The hermiticity relation (4.18) is satisfied by the numerical solutions of (4.8) only up to the 
computing error of the solution of the boundary system (4.9 ). Hence the pencil (4.16 ) is hermitian 
only up to this order. We will force the numerical solution of matrix V to satisfy (4.18) so that 
the hermiticity of the pencil is preserved exactly. This is convenient because the algorithms for 
solving the general eigenvalue problem are much better behaved in the hermitian case |De97| . 

To end this discussion we must realize that, with the basis fa for that we have just con- 
structed, the matrices A and B are almost 3-diagonal and the unique elements different from zero, 
besides the 3-diagonal ones, are those related to the matrix elements of the boundary functions. 
In fact, we can consider a number of cases. If the function fa is an interior bulk function, i.e., 
not corresponding to the nodes X2 nor a;^"^]^, it is obvious that the only nontrivial inner products 
(/a, fb) and {fa, f'b) will correspond to 6 = a — 1, a, a + 1. If the function fa is an extreme bulk 
function, for instance f2°'\ then it has nontrivial inner products only with /^(^a-i) a,nd f^^ ■ If 
fa is now a boundary function fi^^\ then the only non- vanishing inner products will be with the 
other boundary functions and an extreme bulk function, namely /j"'' if / = 2a — 1 or fr°'Li if 
l = 2a. 
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4.5. Convergence of the numerical scheme. We will discuss now the convergence of the 
proposed numerical scheme. The following definitions will help us in this task. 

Definition 4.4. With the previous notations, we will consider the following subspaces: 

i. Du = {^e H\M) I 3f 1^^,^, ^~^^^U{^ + ^^), ^ ^ $|,,,,, ^ = f 

ii. Djj = 1$ G 'H^(M) I if — iip ^ U{(p + iip), = ^l^j^j, 'P = I^I^m}' domain 
of the self-adjoint extension we are interested in. 

Note that Djj C Du is a dense subset of L'^{M) and that the finite dimensional space S*^ — 
span{/a(a;)}, a = 1, . . . , |r|, is contained in Du for all N. 

Lemma 4.5. The closure of Un>qS^ with respect to the Sobolev 1-norm is Djj. 

Proof. The proof uses standard arguments in finite element theory. The family of subdivisions 
of the manifold M defined by the integer N consists on collections ofA'' + n < |r|+n < N + 2n 
closed subintervals Ia,a, o, — 1, . . . , Tq, + 1, of the real line. The domain K of the reference element 
is star-shaped (an interval again) and the family of functions V contains linear polynomials. Thus 
we are in the conditions of Thm. 4.4.4 of |Bre08j and we have the local estimate: 

\\^-^N\\w^^iK) < C(diam(if))™-''^||vl/||,^„(^), 

with 1 < p < oo, and m ~ 1/p > when p > 1. Because we have p — 2, then we consider 
m = 1. The family of subdivisions , consisting on the family of + 1 subintervals la^a of 
length ha of the intervals la, satisfies that max{diam(/Q^a | la. a G M^)} < L/N = h. The 
family of subdivisions is quasi-uniform, i.e., min{diam(/Q^a | la^a G Af^)} > p/{N + 2), if 
< p < mm La. Hence we can apply Thm. 4.4.20 of [Bre08| to our family of subdivisions with 
p = 2, m = 2 and s = 1, to obtain: 

for any 5* G Du, C a constant independent of TV, n the orthogonal projection of to , 
II • ||^2(^,/) the Sobolev 2-norm of ^l* which is finite because Du C H'^{M), and h = L/N as before. 
Then we conclude 

C 

(4.19) 111- - *Ar||«i(M) < J^\\^\hHM)- 

When iV -> cx) we have \\^ - "^NWu^iM) = and therefore = limAr^oo(^'jv) e U'^iM). 
Moreover, because is the linear span of elements in Du, then all functions n satisfy the 
boundary conditions defining Du- As ^ G Du C 'H'^(M) then it has a continuous first derivative 
and a simple use of Sobolev's estimates will show that ^I" G Du too. □ 

In order to prove the converge of the eigenvalues Ajv of the Ritz-Galerkin approximation we 
need a further property of the quadratic form Q. 

Lemma 4.6. The quadratic form Q is continuous with respect to the Sobolev norm \ \ ■ ||i.2- 

Proof. The continuity of the form Q amounts to the continuity of the form a{ip,ip). Notice that 
we can diagonalize the unitary matrix U as: 

U = S^'DS, 

where D = diag(— 1, . . . , —1, zi, . . . , z^), and where Zi are complex numbers of modulus 1 different 
from —1. Now notice that the boundary condition ( |2.6[ ) can be rewritten as 

{I-U)ip^i{I + U)ip. 

Having into account the diagonalization of U we get 

{I - D)Sip = i{I + D)Sip, 
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"2 




V>0 


= i 
























where we have defined Z± = diag(l ± zi, . . . , 1 ± z^) and 



^0 



are just the boundary data adapted to the new basis of ( 



Note that this system of equations 



now imphes that 



3 ipo = and 3 tpi — —iZ_^^Z-ipi. Moreover, as 5 is a unitary 



operator in the Hilbert space of the boundary, the inner product ((^, i^) is given by 
and then 

(4.20) \{^,ip)\<\\Z^Z^'\\\M' <C\ml^. 



It can happen that {—1} ^ o'(J7), then (/ + U) is invertible and equation (4.20) still holds with 
Z± = I ± U. Therefore the quadratic form is continuous with respect to the Sobolev norm 

l|-||l.2. □ 



Now the convergence of the eigenvalues follows by standard arguments. 

Theorem 4.7. The solutions of the approximate eigenvalue problem ($jv,Ajv) converge in the 
limit N ^ oo to the solutions of the weak problem ($, A). 

Proof. The convergence of the eigenvalues is now easily proved by the min-max principle |Re75] , 
which states for the weak eigenvalue problem (4.2 1 that 



sup 

¥'i,...,ip„_i 



inf 



g(<i>,$) 



= A 



and for the approximate eigenvalue problem (4.14) 



sup 

ipi,...,(Pn-l 



inf 



1$ 



N\ 



If we subtract the latter by the former equation we get 

g($Ar,$w) 



(4.21) 



sup 



inf 



N 



sup 



inf 



Q(<i>,<f) 



l$l 



= A 



N 



The l.h.s. goes to zero in the limit — ^ oo by the continuity of the quadratic form Q and therefore 



(4.22) 



lim A^v = A. 

N—^oo 



Finally suppose that $jv is a solution of the approximate eigenvalue problem, i.e. 
(4.23) Qi^N, $jv) = A^v(«'jv, $jv), V*jv G S"" . 

Taking the limit N ^ oo and choosing liniAr^oo '^n — ^ Djj we get 
(4.24) 

hence $ — limjv->.oo $Ar is a solution of the weak problem. 



Q(vl/, lim $jv) = A(5', lim $Ar) V* G Z?;/ 

N-^oo N—¥oo 



□ 
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5. Numerical Experiments 



By now we have proved that the procedure described in section[4]results in the finite dimensional 
eigenvalue problem of the matrix pencil (4.161. Furthermore we have a bound (4.19) on the error 
committed with this discretization. Now if we are able to solve the pencil for increasing grid 
size we will get better and better approximations to the eigenvalue problem. As remarked in 
the previous section, the form of the pencil is almost 3-diagonal and both matrices A and B are 
hermitian matrices, hence the resulting problem is algebraically well behaved. We expose now 
some particular results concerning the stability of the procedure and showing that the bound 
(4.19) is satisfied. We consider for instance the free particle Schrodinger equation in the interval 
[a — 0,b = 2tt], i.e. V{x) = 0. For the sake of simplicity we consider all the physical constants 
normalized to one without loosing generality. In this particular case the matrices A and B take 
the form: 



(5.1a) 



A 



Ml 



-1 



\i21 



il2\ 



2 -1 
-1 A22J 



A- 



2 - Vn -V12 
~V2i 2 - V22 



/-Bu 1/6 
1/6 2/3 



(5.1b) 



B 



Bl2\ 



\B2 



2/3 1/6 
1/6 B22J 



B = 



IV21I 



.^11 



V22V21 + ^12^11] + ^^21 



[V11V12 



22 



V21V22] + 5^12 

+ \Vi2?] + \V2: 



where h 



is the length of each subinterval. Notice that in this simple case (n = 1) there 
is only one interval and we have that r = iV + 1, so we can use r and N almost interchangeably 
if N is large enough, what indeed is the actual case. Now one has to solve for each self-adjoint 
extension the corresponding equations' system (4.8) to obtain the associated pencil. For solving 
the eigenvalue problem of this resulting pencil we have used the Matlab command eig. It uses 
the QZ-algorithm, which is a generalization of the QR-algorithm for the generalized eigenvalue 
problem [Bo07 . In order to check if the bound is satisfied we have to compare the numerical 
solution obtained this way with an analytic one, therefore we will solve the Dirichlet Problem, 
whose solutions are ^{x) — sin ax. In this particular case the solution matrix V is the zero matrix. 
Note that this force the boundary values to be all zero, as it is supposed to be in the Dirichlet 
Problem. We will now observe how 11$ — $7v||«i evolves with increasing N . For calculating it we 
have performed a sum over all the subintervals, knowing that at each subinterval the approximate 
solution is a linear polynomial. In the figure |4] we can see the "H^-norm of the difference between 
the analytic solution for the groundstate of the Dirichlet Problem, i.e. sina;/2, and the computed 
solution plotted against the size of the discretization. 

After a linear fit of the logarithms of this data, one observes that they can be adjusted with 
great goodness to a linear polynomial with slope p = —1.002 ± 0.001, what indeed shows that 



N\\-H- 



= K{r + If < 



K 



1 



< 



K 



therefore satisfying the theoretical bound (4.19) 



Our next step will be to give a proof of the stability of the method against variations of the 
parameters. The unique parameter that this procedure has is the self-adjoint extension whose 
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Figure 4. Evolution of the error, measured in the ?^^-norni, for the groundstate 
of the Dirichlet Problem with increasing lattice size. 



eigenvalue problem we want to solve. We will perturb an initial self-adjoint extension, described 
by the unitary matrix U, and observe the behavior of the relative error in the eigenvalues. In 
other words, we are interested now in studying the relation 

and expect, if the algorithm is stable, that the ratio K(e) grows at most polynomially. However 
one must be very careful in this process since the original eigenvalue problem presents divergences 
in certain circumstances (roughly explained bellow) which could lead to the wrong conclusion 



that this method isn't stable. As a consequence of lemma 4.6 when a self-adjoint extension is 
parameterized by a unitary matrix U that has eigenvalues close but not equal to —1, it may happen 
that the eigenvalues of the considered problem take very large negative values. However, matrices 
with -1 in the spectrum can lead to self-adjoint extensions that are positive definite, for example 
the Dirichlet case. Thus, following a path in the self-adjoint extensions' space, it could happen 
that an infinitesimal change in the arc parameter led to an infinite jump in the eigenvalues. For 
proving stability, it is then necessary to perturb the unitary matrix in a direction of it's tangent 
space in which this jumps don't appear. A path in the self-adjoint extensions' space where these 
jumps don't occur is, for instance, the one described by the so called "quasi-periodic boundary 
conditions" jAs83j . In this case the self-adjoint domain is described by functions that satisfy the 
following boundary conditions u(0) = e*^it(27r) and m'(0) = e'^u'(27r), that in fact correspond to 
the unitary matrix 



m = 



e*^ 

e-^^ 
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Lets then consider perturbations of the periodic case in the quasi-periodic direction, i.e. lets 
consider 

U{e) = U + ieA / 1 ^ ■ ' 1 



1 oy y-i 0^ 

We have calculated the numerical solutions for values of e between 10~^ and 10"'^ in steps of 
10~^ and the discretization size used was N = 250. In this case the relative perturbation is 
||jeA||/||J7|| = e, hence the relative error ratio is 



K{e) 



|A| _ 1 



\\u\\ ' ' 



The results are plotted in figure [sj As can be seen, K{e) grows slowly with e. Fits of the resulting 
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Figure 5. Relative error ratio K{e) of the firsts excited levels for the periodic 
boundary problem plotted against e 



data show that, in every case, the error ratios satisfy power law growths of type K{e) = a • e*" + c 
with exponents -0.89 ± 0.02, -0.42 ± 0.03, -0.03 ± 0.03, 0.28 ± 0.03 for the 1st, 2nd , 3rd and 
4rd excited states respectively. Although the groundstate error ratio has linear growth, and this 
certainly isn't a bad behavior in the sense that it is polynomial, the result is not trustable because 
the quotient |(5A|/|A| is bad conditioned. This is due the fact that the groundlevel in the periodic 
case has A = 0. 



All the numerical calculations of this section have been performed with a laptop computer with 
an Intel Core 2 Duo processor at 2AGHz with 2 Gb RAM. In both examples the calculations 
lasted less than 2 min, which is a further proof of the good discretization procedure, since the eig 
command was called 500 times for the first calculation and 100 times for the second one. 
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